#-----------------------------------------------------------#
# By Keunsoo Kang @ Dankook University (kangk1204@gmail.com)
# Version 1.0.0
# Download listed .sra files and process them
#-----------------------------------------------------------#

use strict;
use warnings;
use Path::Class;

# Directory path to the index
my $SUBREAD_INDEX = "Subread_index/";
my $STAR_INDEX = "STAR_index/";

# General parameters
my $CPU = 4;	# Number of threads

# Files
my @GSM_INFO;	# .octo list file

# Commands
my $SRA_CMD = "fastq-dump --split-3 ";	# SRATool kit command: SRA to FASTQ
my $ASPERA_CMD = "~/.aspera/connect/bin/./ascp -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh -k 1 -Tr -l500m anonftp\@ftp-private.ncbi.nlm.nih.gov:";	# SRA files download by Aspera
my $PHRED_CMD = " \| head \-n1000000 \| awk \'\{if\(NR\%4\=\=0\) printf\(\"\%s\"\,\$0\)\;\}\' \|  od \-A n \-t u1 \| awk \'BEGIN\{min\=100\;max\=0\;\}\{for\(i\=1\;i\<\=NF\;i\+\+\) \{if\(\$i\>max\) max\=\$i\; if\(\$i\<min\) min\=\$i\;\}\}END\{if\(max\<\=74 \&\& min\<59\) print \"Phred\+33\"\; else if\(max\>73 \&\& min\>\=64\) print \"Phred\+64\"\; else if\(min\>\=59 \&\& min\<64 \&\& max\>73\) print \"Solexa\+64\"\; else print \"Unknown score encoding\!\"\;\}\' \> ";	# Phred check; Thanks to Xianjun Dong (http://onetipperday.blogspot.kr/2012/10/code-snip-to-decide-phred-encoding-of.html)
my $SUBR_CMD = "subread-align -u -i $SUBREAD_INDEX";	# Mapping command by Subread
my $STAR_CMD = "STAR --outSAMstrandField intronMotif --runThreadN $CPU --genomeDir $STAR_INDEX";	# RNA-seq mapping by STAR
my $CHR_SIZE_CMD = "fetchChromSizes ";
my $HOMERDIR_CMD = "makeTagDirectory ";	# makeTagDirectory by HOMER
my $HOMERRM_CMD = "removeOutOfBoundsReads.pl ";	# Remove out range reads by HOMER
my $HOMERBW_CMD = "makeUCSCfile ";	# makeUCSCfile bigwig by HOMER

# Directories predefined by the update_the_list.pl script

my $SRR_TMP_DIR = "00_SRA_files"."\/";
my $HOMER_DIR = "03_Subread_HOMER_dir\/";
my $INDEL_DIR = "04_Subread_INDEL_dir\/";



#-------------------------------------[ Main started ]-------------------------------------#
GET_FILE_INFO();
DOWNLOAD_N_PROCESS();




#---------------------------------------[ Functions ]--------------------------------------#
sub GET_FILE_INFO{
	# Grab the .octo list file from the 02_LIST_dir directory
	dir('.')->recurse(callback => sub{ 
		my $file = shift;
		if($file =~ /\.octo$/){
			push my @files, $file->absolute->stringify;

			for my $f(0..$#files){
				if($files[$f] =~ /02_LIST_dir/gi){
					open(FO, $files[$f]) or die $!;
					while(<FO>){
						chomp;

						push( @GSM_INFO, $_ );
					}
					close(FO);
				}
			}
		}
	});		
}

sub DOWNLOAD_N_PROCESS{
	for my $f(0..$#GSM_INFO){
		my @d = split/\t/, $GSM_INFO[$f];
		my $gse = $d[1];
		my $gsm = $d[2];
		my $pmid = $d[7];
		my $type = $d[5];	# ChIP-seq/MNase-seq or RNA-seq?
		my $index = $d[4];	# Genome assembly?

		my $header = "$index\_$gsm\_$gse\_$pmid\_$type";

		if($index =~ /unknown/gi){
			timestamp();
			print "# Processing.. $header\n";
			timestamp();
			print "   Unable to process further due to invalid genome index: $index\n";
			next;
		}

		# FILTER by type and genome assembly
		if($type =~ /ChIP-seq|MNase-seq|DNase-seq/gi){
			generate_subread_index($index);

			timestamp();
			print "# Processing.. $header\n";
			non_rna_seq($GSM_INFO[$f], $header);	# By HOMER
		}elsif($type =~ /RNA-seq/gi){	
			generate_star_index($index);

			timestamp();
			print "# Processing.. $header\n";
			rna_seq($GSM_INFO[$f], $header);	# By STAR
		}else{
			timestamp();
			print "# Processing.. $header\n";
			timestamp();
			print "   Unable to process further due to invalid NGS data type: $type -> SKIPPED!\n";
			next;
		}
	}
}

sub generate_subread_index{
	my $genome = shift;

	if(-e $SUBREAD_INDEX.$genome."\.00\.b\.tab"){
	}else{
		unless(-d $SUBREAD_INDEX){
			my $mkdir_cmd = "mkdir $SUBREAD_INDEX";
			system $mkdir_cmd;
		}

		unless(-e $genome."\.fa"){
			my $wget_cmd = "wget http://143.248.14.46/Octopus/".$genome."\.fa";
			timestamp();
			print "! Downloading $genome sequence from the Octopus server..\n";
			system $wget_cmd;
		}

		my $subread_cmd = "subread-buildindex -o ".$SUBREAD_INDEX.$genome." $genome\.fa";
		timestamp();
		print "! Generating subread index for $genome in $SUBREAD_INDEX\n";
		system $subread_cmd;
	}
}

sub generate_star_index{
	my $genome = shift;

	if(-e $STAR_INDEX.$genome."\/SAindex"){
	}else{
		unless(-d $STAR_INDEX.$genome){
			my $mk_dir = "mkdir -p $STAR_INDEX".$genome;
			system $mk_dir;
		}

		unless(-e $genome."\.fa"){
			my $wget_cmd = "wget http://143.248.14.46/Octopus/".$genome."\.fa";
			timestamp();
			print "! Downloading $genome sequence from the Octopus server..\n";
			system $wget_cmd;
		}

		my $star_cmd = "STAR --runMode genomeGenerate --runThreadN $CPU --genomeDir ".$STAR_INDEX.$genome." --genomeFastaFiles $genome\.fa";
		timestamp();
		print "! Generating star index for $genome in $STAR_INDEX\n";
		system $star_cmd;
	}
}

sub rna_seq{
	my $detail = shift;
	my $header = shift;

	my @d = split/\t/, $detail;	

	my $gse = $d[1];
	my $gsm = $d[2];
	my $srp = $d[3];
	my $index = $d[4];
	my $type = $d[5];	# RNA-seq?
	my $single_pair = $d[6];
	my $pmid = $d[7];
	my $desc = $d[10];
	my @sra_file = split/\,/, $d[8];
	my @ftp_link = split/\,/, $d[9];

	my @sam_files;
	my $chrom_size_file = $index."\.chrom.sizes";

	# Check the genome index file
	unless(-e $STAR_INDEX.$index."\/SAindex"){
		timestamp();
		print "  No genome index file for STAR -> ".$STAR_INDEX.$index."\/SAindex etc..SKIPPED!\n";
		return;
	}	

	# Data download
	for my $s(0..$#sra_file){
		my @tt = split/ftp\:\/\/ftp\-trace\.ncbi\.nlm\.nih\.gov\//, $ftp_link[$s];
		my $aspera_cmd = $ASPERA_CMD.$tt[1].$sra_file[$s]."\/".$sra_file[$s]."\.sra ".$SRR_TMP_DIR.$sra_file[$s]."\.sra";
		timestamp();
		print "  Downloading.. ".$sra_file[$s]."\.sra by aspera\n";
		system $aspera_cmd;	
	}

	# Data processing
	for my $s(0..$#sra_file){

		my $single_or_paired = 0;

		my $files = $SRR_TMP_DIR.$sra_file[$s]."\.sra";

		my @tmpone = split/\//, $files;
		my @tmptwo = split/\./, $tmpone[$#tmpone];
		my $filealias = $tmptwo[0];

		# SRATool kit
		timestamp();
		print "  Converting $files to fastq by SRATool kit\n";
		my $sra_to_fastq_cmd = $SRA_CMD.$files;
		system $sra_to_fastq_cmd;

		my $phred_status;

		# Determine single- or paired-end strategy
		if(-e $filealias."\_1\.fastq"){
			$single_or_paired = 2;

			# Determine phred scheme
			timestamp();
			$phred_status = which_phred($filealias."\_1\.fastq");
			print "  Determine the phred scheme -> $phred_status\n";
		}elsif(-e $filealias."\.fastq"){
			$single_or_paired = 1;

			# Determine phred scheme
			timestamp();
			$phred_status = which_phred($filealias."\.fastq");
			print "  Determine the phred scheme -> $phred_status\n";
		}else{
			$single_or_paired = 0;
			log_stamp($files, "\.sra2fastq_failure\.octolog");
			timestamp();
			print "! Conversion failure... please check the $files file manually\n";
			next;
		}						

		# STAR
		timestamp();
		print "  Mapping $filealias fastq file(s) to $index genome by STAR\n";
		my $subread_cmd;
		my $phred_option;
		if($phred_status =~ /Phred\+33/){
			$phred_option = "-P 3";
		}elsif($phred_status =~ /Phred\+64/){
			$phred_option = "-P 6";
		}elsif($phred_status =~ /Solexa\+64/){
			$phred_option = "na";
			# Skipping
		}else{
			$phred_option = "na";
			# Skipping
		}

		if($single_or_paired == 2){
			if($phred_option eq "na"){
				timestamp();
				print "! Cannot detect the phred score.. skipped\n";

				my $rm_fastq_cmd = "rm $filealias\_1\.fastq";
				system $rm_fastq_cmd;
				$rm_fastq_cmd = "rm $filealias\_2\.fastq";
				system $rm_fastq_cmd;

				last;
			}

#			my $subread_cmd = $SUBR_CMD.$index."sub ".$phred_option." -r ".$filealias."\_1\.fastq"." -R ".$filealias."\_2\.fastq -o pe\_$filealias\.sam";
			my $star_cmd = $STAR_CMD.$index." --readFilesIn ".$filealias."\_1\.fastq $filealias"."\_2\.fastq";
			timestamp();
			print "  Paired-end mapping started\n";
			system $star_cmd;
			log_stamp($files, "\.STARmapping_completed\.octolog");

			# Change the aligned file namd
			my $change_cmd = "mv Aligned.out.sam pe\_$filealias\.sam";
			system $change_cmd;

			push( @sam_files, "pe\_$filealias\.sam" );

			my $rm_fastq_cmd = "rm $filealias\_1\.fastq";
			system $rm_fastq_cmd;
			$rm_fastq_cmd = "rm $filealias\_2\.fastq";
			system $rm_fastq_cmd;						

		}elsif($single_or_paired == 1){
			if($phred_option eq "na"){
				timestamp();
				print "! Cannot detect the phred score.. skipped\n";

				my $rm_fastq_cmd = "rm $filealias\.fastq";
				system $rm_fastq_cmd;

				last;
			}

#			my $subread_cmd = $SUBR_CMD.$index."sub ".$phred_option." -r ".$filealias."\.fastq"." -o se\_$filealias\.sam";
			my $star_cmd = $STAR_CMD.$index." --readFilesIn ".$filealias."\.fastq";
			timestamp();
			print "  Single-end mapping started\n";
			system $star_cmd;
			log_stamp($files, "\.STARmapping_completed\.octolog");

			# Change the aligned file namd
			my $change_cmd = "mv Aligned.out.sam se\_$filealias\.sam";
			system $change_cmd;

			push( @sam_files, "se\_$filealias\.sam" );

			my $rm_fastq_cmd = "rm $filealias\.fastq";
			system $rm_fastq_cmd;
		}
	}

	if(@sam_files == @sra_file){	# All the sra files were mapped successfully,
		# Delete the sra files and generate a code file
		for my $f(0..$#sra_file){
			my $rm_cmd = "rm ".$SRR_TMP_DIR.$sra_file[$f]."\.sra";
			system $rm_cmd;
		}

		# Chromosome size file
		unless(-e $index."\.chrom.sizes"){
			timestamp();
			my $fetch_cmd = $CHR_SIZE_CMD."$index \> $index"."\.chrom\.sizes";
			print "  Generating chromosome size file -> $index\.chrom\.sizes";
			$chrom_size_file = $index."\.chrom\.sizes";
			system $fetch_cmd;
		}

		# HOMER analysis
		my $homer_file_input = join(" ", @sam_files);

		timestamp();
		my $tagdir_cmd = $HOMERDIR_CMD.$HOMER_DIR."$index\/$header\/ ".$homer_file_input." -format sam -fragLength given";
		print "# SAMtoBigWig: $tagdir_cmd\n";
		system $tagdir_cmd;

		my $remove_cmd = $HOMERRM_CMD.$HOMER_DIR."$index\/$header\/ none -chromSizes $chrom_size_file";
		system $remove_cmd;

		my $bigwig_cmd = $HOMERBW_CMD.$HOMER_DIR."$index\/$header\/ -fragLength given -fsize 5e7 -bigWig $chrom_size_file -o auto >> trackinfo.txt";
		system $bigwig_cmd;
	}

	# Remove the sam file
	for my $k(0..$#sam_files){
		my $rmsam_cmd = "rm $sam_files[$k]";
		system $rmsam_cmd;
	}
}

sub non_rna_seq{
	my $detail = shift;
	my $header = shift;

	my @d = split/\t/, $detail;	

#	print "@d\n";

	my $gse = $d[1];
	my $gsm = $d[2];
	my $srp = $d[3];
	my $index = $d[4];
	my $type = $d[5];	# ChIP-seq?
	my $single_pair = $d[6];
	my $pmid = $d[7];
	my $desc = $d[10];
	my @sra_file = split/\,/, $d[8];
	my @ftp_link = split/\,/, $d[9];

	my @sam_files;
	my $chrom_size_file = $index."\.chrom.sizes";

	# Check the genome index file
	unless(-e $SUBREAD_INDEX.$index."\.00\.b\.tab"){
		timestamp();
		print "  No genome index file for Subread -> $SUBREAD_INDEX".$index."\.00\.b\.tab etc..SKIPPED!\n";
		return;
	}
	
	# Data download
	for my $s(0..$#sra_file){
		my @tt = split/ftp\:\/\/ftp\-trace\.ncbi\.nlm\.nih\.gov\//, $ftp_link[$s];
		my $aspera_cmd = $ASPERA_CMD.$tt[1].$sra_file[$s]."\/".$sra_file[$s]."\.sra ".$SRR_TMP_DIR.$sra_file[$s]."\.sra";
		timestamp();
		print "  Downloading.. ".$sra_file[$s]."\.sra by aspera\n";
		system $aspera_cmd;	
	}

	# processing
	for my $s(0..$#sra_file){

		my $single_or_paired = 0;

		my $files = $SRR_TMP_DIR.$sra_file[$s]."\.sra";

		my @tmpone = split/\//, $files;
		my @tmptwo = split/\./, $tmpone[$#tmpone];
		my $filealias = $tmptwo[0];

		# SRATool kit
		timestamp();
		print "  Converting $files to fastq by SRATool kit\n";
		my $sra_to_fastq_cmd = $SRA_CMD.$files;
		system $sra_to_fastq_cmd;

		my $phred_status;

		# Determine single- or paired-end strategy
		if(-e $filealias."\_1\.fastq"){
			$single_or_paired = 2;

			# Determine phred scheme
			timestamp();
			$phred_status = which_phred($filealias."\_1\.fastq");
			print "  Determine the phred scheme -> $phred_status\n";
		}elsif(-e $filealias."\.fastq"){
			$single_or_paired = 1;

			# Determine phred scheme
			timestamp();
			$phred_status = which_phred($filealias."\.fastq");
			print "  Determine the phred scheme -> $phred_status\n";
		}else{
			$single_or_paired = 0;
			log_stamp($files, "\.sra2fastq_failure\.octolog");
			timestamp();
			print "! Conversion failure... please check the $files file manually\n";
			next;
		}						

		# Subread
		timestamp();
		print "  Mapping $filealias fastq file(s) to $index genome by Subread\n";
		my $subread_cmd;
		my $phred_option;
		if($phred_status =~ /Phred\+33/){
			$phred_option = "-P 3";
		}elsif($phred_status =~ /Phred\+64/){
			$phred_option = "-P 6";
		}elsif($phred_status =~ /Solexa\+64/){
			$phred_option = "na";
			# Skipping
		}else{
			$phred_option = "na";
			# Skipping
		}

		if($single_or_paired == 2){
			if($phred_option eq "na"){
				timestamp();
				print "! Cannot detect the phred score.. skipped\n";

				my $rm_fastq_cmd = "rm $filealias\_1\.fastq";
				system $rm_fastq_cmd;
				$rm_fastq_cmd = "rm $filealias\_2\.fastq";
				system $rm_fastq_cmd;

				last;
			}

			my $subread_cmd = $SUBR_CMD.$index." ".$phred_option." -r ".$filealias."\_1\.fastq"." -R ".$filealias."\_2\.fastq -o pe\_$filealias\.sam";
			timestamp();
			print "  Paired-end mapping started\n";
			system $subread_cmd;
			log_stamp($files, "\.mapping_completed\.octolog");
			push( @sam_files, "pe\_$filealias\.sam" );

			my $rm_fastq_cmd = "rm $filealias\_1\.fastq";
			system $rm_fastq_cmd;
			$rm_fastq_cmd = "rm $filealias\_2\.fastq";
			system $rm_fastq_cmd;						

			# INDEL file
			my $mv_indel = "mv pe\_$filealias\.sam\.indel $INDEL_DIR";
			system $mv_indel;

		}elsif($single_or_paired == 1){
			if($phred_option eq "na"){
				timestamp();
				print "! Cannot detect the phred score.. skipped\n";

				my $rm_fastq_cmd = "rm $filealias\.fastq";
				system $rm_fastq_cmd;

				last;
			}

			my $subread_cmd = $SUBR_CMD.$index." ".$phred_option." -r ".$filealias."\.fastq"." -o se\_$filealias\.sam";
			timestamp();
			print "  Single-end mapping started\n";
			system $subread_cmd;
			log_stamp($files, "\.mapping_completed\.octolog");
			push( @sam_files, "se\_$filealias\.sam" );

			my $rm_fastq_cmd = "rm $filealias\.fastq";
			system $rm_fastq_cmd;

			# INDEL file
			my $mv_indel = "mv se\_$filealias\.sam\.indel $INDEL_DIR";
			system $mv_indel;
		}
	}

	if(@sam_files == @sra_file){	# All the sra files were mapped successfully,
		# Delete the sra files and generate a code file
		for my $f(0..$#sra_file){
			my $rm_cmd = "rm ".$SRR_TMP_DIR.$sra_file[$f]."\.sra";
			system $rm_cmd;
		}

		# Chromosome size file
		unless(-e $index."\.chrom.sizes"){
			timestamp();
			my $fetch_cmd = $CHR_SIZE_CMD."$index \> $index"."\.chrom\.sizes";
			print "  Generating chromosome size file -> $index\.chrom\.sizes";
			$chrom_size_file = $index."\.chrom\.sizes";
			system $fetch_cmd;
		}

		# HOMER analysis
		my $homer_file_input = join(" ", @sam_files);

		timestamp();
		my $tagdir_cmd = $HOMERDIR_CMD.$HOMER_DIR."$index\/$header\/ ".$homer_file_input." -format sam";
		print "# SAMtoBigWig: $tagdir_cmd\n";
		system $tagdir_cmd;

		my $remove_cmd = $HOMERRM_CMD.$HOMER_DIR."$index\/$header\/ none -chromSizes $chrom_size_file";
		system $remove_cmd;

		my $bigwig_cmd = $HOMERBW_CMD.$HOMER_DIR."$index\/$header\/ -fsize 5e7 -bigWig $chrom_size_file -o auto >> trackinfo.txt";
		system $bigwig_cmd;
	}

	# Remove the sam file
	for my $k(0..$#sam_files){
		my $rmsam_cmd = "rm $sam_files[$k]";
		system $rmsam_cmd;
	}
}

sub log_stamp{
	my $path = shift;
	my $output = shift;

	my @tmp = split/\//, $path;
	my $loc = join("\/", @tmp[0..$#tmp]);
	open(FFOUT, ">".$loc.$output) or die $!;
	close(FFOUT);	
}

sub which_phred{
	my $fastq_file = shift;

	my $rand_num = int(rand(100000));
	my $ph_cmd = "less $fastq_file".$PHRED_CMD."$rand_num\.txt";
	system $ph_cmd;

	my $phred_type;
	open(FRR, "$rand_num\.txt") or die $!;
	while(<FRR>){
		chomp;
		$phred_type = $_;
		last;
	}
	close(FRR);	

	my $rm_command = "rm $rand_num\.txt";
	system $rm_command;

	return $phred_type;
}




#---------------------------------------------------------------[ MISC functions ]
sub timestamp{
	my ($sec, $min, $hour, $mday, $mon, $year, $wday, $yday, $isdst) = localtime(time);
	printf "%4d-%02d-%02d %02d:%02d:%02d ", $year+1900,$mon+1,$mday,$hour,$min,$sec;
}

#	# Grab the SRA files
#	dir('.')->recurse(callback => sub{ 
#		my $file = shift;
#		if($file =~ /\.sra/){
#			push @SRA_FILES, $file->absolute->stringify;
#		}
#	});

sub uniqueElements {
	my ($item, %seen, @result);
	foreach $item (@_) {
		push(@result, $item) unless $seen{$item}++;
	}
	
	return @result;
}
